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ABSTRACT 

Granular  Resistive  Force  Theory  (RFT)  is  a  reduced-order  model  inspired  by  analogous  boundary-integral  methods  for  Stokesian  fluids. 
Despite  its  remarkable  capability  to  predict  experimental  locomotion  and  force  distributions  on  mobile  bodies  in  granular  media,  there  is  no 
theoretical  understanding  for  this  behavior.  Moreover,  such  a  reduction  is  surprising  given  the  highly  nonlinear  constitutive  behavior  of 
granular  media.  There  could  be  a  variety  of  reduced-order  applications  in  various  structure/granular  interactions  if  a  theoretical  picture  in 
explanation  of  RFT  could  be  uncovered.  Moreover,  as  with  any  accurate  and  sufficiently  reduced  model,  RFT,  once  its  limitations  and 
backing  is  better  understood,  could  be  used  as  a  reliable  design-optimization  tool  to  produce  locomotors  with  optimal  shapes,  or  tunable 
flexibility  to  improve  efficiency  of  locomotion  within  granular  media. 

We  have  performed  initial  investigations  of  these  questions.  A  systematic  set  of  continuum  simulations  were  conducted.  In  particular,  we 
have  determined  that  the  RFT  superposition  principle  arises  naturally  from  full-field  solutions  of  plasticity  under  Drucker-Prager  yielding, 
the  most  commonly  used  plastic  flow  model  for  granular  media  in  engineering.  We  find  that  the  empirical  input  data  used  in  RFT,  fit  from 
many  experiments,  in  fact  arises  quantitatively  from  plastic  flow  solution,  requiring  only  a  friction  coefficient  and  material  density  to  solve. 
The  results  have  motivated  mathematical  analyses  to  identify  the  approximations  at  play  to  achieve  the  desired  model  reduction.  In  parallel, 
a  design-optimization  protocol  for  wheeled  locomotors  in  granular  beds  has  been  studied,  based  on  numerical  optimization  of  wheel  shape. 
This  effort  has  shown  that  under  certain  conditions,  the  natural  curvature  of  a  wheel  tread  under  axle-bearing  weight  is  proper  to  improve 
traction  in  soils. 
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Abstract 

The  flowability  of  loose  terrain  under  solid  intrusion  produces  complex  dynamics  in  problems 
ranging  from  geotechnical  design  to  animal  and  vehicle  locomotion.  One  approach  is  the  Resistive 
Force  Theory  (RFT),  a  recent  empirical  tool  for  predicting  granular-structure  interaction.  Its 
simplicity  and  effectiveness  are  surprising  given  the  many  complexities  of  granular  flow,  begging 
fundamental  questions  of  why  RFT  works.  We  have  found  a  link  between  RFT  and  plasticity 
theory,  showing  RFT  arises  solely  from  frictional  yielding.  Without  any  fitting,  plasticity  generates 
experimental  RFT  data,  and  reproduces  RFT’s  foundational  assumptions  including  spatial  force 
superposition.  Continuum  dimensional  analysis  explains  why  RFT  is  more  accurate  in  granular 
media  than  its  older  counterpart  for  viscous  fluids,  which  instructs  the  development  of  RFT’s  in 
other  media  and  for  more  diverse  conditions. 


One  Sentence  Summary: 

We  have  discovered  a  fundamental  link  between  the  theory  of  granular  plasticity  and  the  as-yet 
empirical  Resistive  Force  Theory  of  terradynamics. 

Main  Text: 

Introduction  -The  interaction  of  solid  objects  with  granular  media  is  a  common  aspect  of  many 
engineering  processes  as  seen  in  several  examples  of  bio-locomotion  and  hydrodynamics  [1,  2].  Inter¬ 
estingly,  such  a  broad  problem  is  not  trivial  to  solve  in  granular  materials  due  to  the  complexity  of 
the  response  of  the  system  and  the  many  interaction  forces  that  define  the  movement  at  the  grain 
level,  leaving  explorations  in  this  area  to  rely  upon  experimentally  driven  studies  or  numerical  meth¬ 
ods  such  as  Discrete  Element  Method  (DEM).  Continuum  modeling  techniques  that  can  accurately 
represent  collective  grain  interactions  on  the  macroscopic  scale  would  be  more  efficient  than  DEM 
computationally,  providing  a  route  toward  optimization  in  granular-structure  interaction,  and  could 
yield  new  fundamental  insights  into  the  physics  of  granular  resistance. 

On  the  topic  of  interaction  of  granular  material  with  intruding  objects,  attempts  have  been  made 
mostly  driven  by  experimental  observations  to  describe  this  interaction.  Despite  the  fact  that  a 
fundamental  derivation  is  missing,  a  simple  yet  very  effective  empirical  tool  known  as  Resistive  Force 
Theory  (RFT)  for  granular  materials  has  been  proposed  to  approximate  the  forces  on  intruding 
objects  moving  through  granular  media.  When  coupled  to  the  momentum  balance  equations,  it 
provides  a  simple  and  predictive  tool  for  simulating  the  locomotion  of  arbitrarily  shaped  moving 
bodies  in  loose  terrain  [2-4].  The  simplicity  of  the  theory  and  its  predictive  effectiveness  are  surprising 
in  light  of  the  complex,  nonlinear,  and  oftentimes  visibly  nonlocal  constitutive  properties  of  granular 
media  [5-9]. 
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A  fundamental  problem,  then,  is  to  determine  why  RFT  works.  Beyond  theoretical  concerns, 
understanding  the  origin  of  this  rule  could  shed  light  on  the  requirements  for  its  validity  and  the  limits 
of  its  usage.  In  this  paper  we  provide  an  answer  to  this  question  by  supplying  numerical  and  analytical 
evidence  that  RFT  arises  from  the  partial  differential  equations  (PDE’s)  of  frictional  plasticity  theory. 
To  demonstrate  this  connection,  we  focus  on  dry  granular  material  which  is  compliant  to  laboratory 
experiment  and  is  representative  of  a  wide  class  of  flowable  materials.  Resistive  Force  Theory  -  RFT 
was  initially  developed  to  approximate  the  speed  of  swimming  micro-organisms  at  low  Reynolds 
numbers  [10]  by  studying  the  thrust  and  drag  of  individually  moving  elements  of  its  body.  This 
method  has  been  used  extensively  to  study  the  hydrodynamics  and  mobility  of  cells  as  reviewed  in 
[11].  In  the  study  of  granular  media,  a  similar  notion  arose  much  more  recently,  when  in  experimental 
studies  of  arbitrarily-shaped  intruders  moving  in  granular  beds,  it  was  determined  that  the  resistive 
force  against  intruder  motion  is  rather  well  represented  by  a  simple  superposition  rule  [4];  the  intruder 
boundary  can  be  decomposed  into  a  connected  collection  of  differential  planar  elements  and  the  total 
resistive  force  is  deemed  equal  to  the  sum  of  the  resistive  forces  on  each  plane  as  if  it  was  moving 
steadily  on  its  own.  RFT  can  be  used  in  general  3D  setups,  but  for  concreteness,  let  us  consider  now 
the  case  of  a  quasi-2D  intruder  having  some  thickness  D  which  displaces  in  the  xz- plane  (gravity 
points  in  —  £  direction,  z  =  0  represents  the  sand  surface)  and  whose  cross-section  also  lies  in  the 
xz- plane.  For  any  subset  S  of  the  leading  surface  of  the  intruder,  RFT  states  the  resistive  force  on 
S  will  obey 

(. fx,fz)=  [  (ax(/3, 'j),  az(l3, 7))  H(z)  \z\  dS  (1) 

Js 

where  /?  is  the  orientation  angle  (attack  angle)  of  the  differential  surface  element  and  7  is  the  angle  of 
the  velocity  angle  (intrusion  angle)  of  the  surface  element,  both  measured  from  the  horizontal.  The 
Heaviside  function,  H(z),  removes  resistive  force  above  the  granular  level  and  increases  resistance 
with  depth  to  reflect  the  increase  in  granular  pressure  with  depth.  The  key  constitutive  ingredient 
in  the  theory  is  the  selection  of  the  two  functions  otx  and  07,  which  is  done  empirically  by  fitting 
experimental  force  data  on  intruding  flat  plates  under  various  7  and  f3  conditions. 

Continuum  Modeling  -  Frictional  plasticity  models  based  on  the  Mohr-Coulomb  (MC)  yield  cri¬ 
terion  or  Drucker-Prager  (DP)  criterion  are  the  most  commonly  used  granular  continuum  models  for 
predicting  the  deformation  of  a  body  of  grains.  We  choose  to  study  the  DP  yield  criterion  in  this 
paper  as  it  is  simpler  to  implement  in  3D,  and  in  2D  the  theory  matches  MC.  We  will  make  two 
important  choices  in  closing  the  system  of  constitutive  equations:  (i)  we  choose  a  non-associative 
flow  rule  corresponding  to  incompressible  flow  during  yielding,  to  emulate  an  immediate  approach 
to  a  critical  state  [5]  of  volume-conserving  flow,  and  (ii)  we  append  an  ‘opening  rule’  to  the  system 
to  remove  granular  stress  when  a  material  element  attempts  to  enter  a  state  of  tension,  as  granular 
media  naturally  disconnects  when  a  tensile  state  is  attempted.  Open  states  often  occur  temporarily 
in  the  wake  of  a  moving  intruder,  before  material  from  above  collapses  onto  the  material  back  down 
(see  Supplemental  Materials  for  more  details).  The  bare-bones  granular  flow  model  described  above 
does  not  account  for  rate-sensitivity,  strain-dependent  strength/porosity,  fabric  anisotropy  during 
flow,  or  nonlocal  effects  based  on  particle  size,  which  are  all  known  to  exist. 

The  strain-rate  tensor  is  defined  from  the  spatial  velocity  field,  77,  by  D\j  —  (dvi/ dxj+dvj / dxi) / 2. 
We  define  the  scalar  (equivalent)  shear-rate  as  7  =  ^j2D'i-Dri-  where  D[-  —  Dij  —  SijDu/3  is  the 
strain-rate  deviator.  Assuming  that  the  Cauchy  stress,  07/ ,  is  co-directional  with  the  strain-rate,  the 
DP  criterion  reduces  to 

Gij  =  -PSij  +  2/jLcPDfij/'y  if  7,  P  >  0.  (2) 

In  the  above,  fic  is  a  constant  internal  friction  coefficient  and  P  is  the  isotropic  pressure.  Whenever 


10  cm/s 


Figure  1:  Speed  (left  two  columns)  and  velocity  directions  (right  two  columns)  due  to  the  motion  of  a  buried  flat 
intruder  of  different  orientations  moving  to  the  right.  Top  row:  FEM  simulations  of  the  continuum  model.  Bottom 
row:  DEM  simulations  of  [3].  Velocity  arrows  corresponding  to  regions  moving  at  the  speed  of  more  than  5  cm/s  are 
as  red. 


7,  P  >  0  we  assume  incompressible  plastic  flow  (Du  =  0)  such  that  the  material  density  p  remains 
at  a  chosen  ‘dense’  value  pc.  Whenever  p  <  pc,  we  set  oy/  =  0  to  represent  granular  disconnection. 
To  provide  a  stress  in  rigid  zones,  where  7  =  0  and  P  >  0,  and  to  aid  in  implementing  the  pressure 
constraints  numerically,  we  include  a  small  elastic  part  to  the  deformation  [12];  as  long  as  the  elastic 
stiffness  is  sufficiently  high  the  observed  plastic  flow  behavior  is  unaffected,  a  point  we  verified 
directly  in  our  simulations.  Momentum  balance,  daij/dxj  +  pgi  =  piy  closes  the  system  for  arbitrary 
boundary  value  problems,  where  gi  is  the  acceleration  of  gravity  and  the  superscript  dot  represents 
the  Lagrangian  time  derivative. 

FEM  Results  -  The  model  just  described  is  numerically  implemented  in  3D  using  Abaqus/ Explicit 
[13].  The  constitutive  relation  is  inputted  as  a  custom  VUMAT  subroutine.  We  first  consider 
problems  with  plane-strain  symmetry.  We  use  second  order  explicit  tetrahedral  elements  in  Abaqus 
package  to  discretize  the  granular  material  space  into  a  finite  element  mesh.  The  meshed  space  is  a 
row  of  single  element  thickness  and  consequently  resembles  plane  strain  condition.  On  any  container 
edge,  no-penetration  conditions  are  applied.  Gravity  is  applied  to  obtain  the  proper  depth-dependent 
pressure  distribution  before  intruder  motion  begins.  We  model  the  intruder  as  a  rigid  surface  within 
the  granular  body  by  assigning  corresponding  planar  nodes  to  move  as  a  rigid  body  at  a  constant 
rate.  This  resembles  a  fully  rough  surface  due  to  the  no-slip  condition,  no  penetration  condition.  We 
use  sufficiently  refined  second-order  elements  to  maintain  solution  accuracy.  Sample  flow  directions 
and  velocity  distributions  within  the  granular  material  obtained  by  FEM  are  shown  in  Fig.  1  and 
compared  to  the  DEM  results  in  the  literature  [3].  It  is  worthwhile  to  note  that  at  the  intruder 
interface,  grains  naturally  tend  to  slide  or  rotate,  so  our  no-slip  condition  is  hard  to  physically 
replicate.  But  this  approach  relinquishes  the  need  to  define  a  complex  contact  routine  which  brings 
in  additional  contact  parameters. 

To  perform  a  quantitative  analysis  of  the  approach  and  to  establish  the  connection  between 
frictional  plasticity  and  RFT,  we  repeated  the  method  described  above  on  an  intruder  at  attack 
angles  (5  and  intrusion  angles  7  (0  <  /3,y  <  7r).  The  drag  and  lift  forces  acting  on  the  plate  are 


Figure  2:  RFP’s  obtained  from  FEM  simulations  with  frictional  plasticity  (‘sim.’  superscript),  compared  against 
experimental  RFP’s  (‘exp.’  superscript)  from  [2]  .  Solid  black  line  shows  the  zero  contour. 


then  converted  to  the  two  functions  ax  and  az,  which  is  done  by  averaging  over  a  period  of  time 
after  plastic  deformation  is  well  developed.  We  use  the  properties  of  loose  packed  3mm  glass  beads 
as  reported  in  [2]:  the  friction  coefficient  is  /ic  =  0.38  based  on  their  repose  angle  calculation,  the 
density  is  p  =  2.6  g/cm3,  and  the  packing  fraction  of  0.63.  As  plotted  in  Fig. 2,  the  resistive  force 
plots  (RFP)  of  olx  and  az  are  strikingly  similar  to  the  experimentally  obtained  RFPs  reported  in  [2] 
for  loose-packed  3mm  glass  beads  that  show  quantitative  agreement.  Furthermore  the  zero  contour 
lines  in  the  numerical  RFPs  follow  the  same  path  as  that  of  the  experimental  RFPs.  The  only 
noticeable  difference  between  the  two  sets  of  figures  is  observed  in  the  location  of  maximum  drag 
force  in  the  )  plot,  which  can  be  attributed  to  the  no-slip  boundary  condition  used  in  the 

simulations.  What  is  more  remarkable  about  this  figure  is  that  no  fitting  parameters  are  used  in  the 
constitutive  model,  only  the  physically  measurable  quantities  of  density,  repose  angle  and  packing 
fraction. 

These  plots  are  also  noticeably  accurate  in  reproducing  the  resistive  force  distribution  of  an 
arbitrarily  shaped  object  moving  within  continuum  granular  material.  For  instance,  in  Fig.  3,  the 
distribution  and  resultant  forces  acting  on  circular  and  diamond-shaped  intruders  is  obtained  directly 
by  FEM  and  compared  with  those  predicted  by  the  FEM-generated  RFP’s  under  the  form  of  Eq. 
1.  These  illustrations  show  the  extent  of  the  validity  of  the  superposition  rule  of  RFT  and  further 
show  the  connection  between  frictional  plasticity  and  RFT.  Though  some  errors  are  observed  at  the 
edges,  the  distribution  and  total  forces  are  closely  comparable.  Normal  and  shear  stresses  on  the 
perimeter  of  the  leading  edge  of  the  circular  object  show  a  very  close  match  in  the  bottom  quarter, 
but  deviations  pick  up  in  the  top  quarter  which  is  attributed  to  the  no-slip  condition.  The  material  at 
the  backward  inclined  surface  tends  to  slide  in  reality,  where  in  case  of  our  simulation  this  movement 
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Figure  3:  Distribution  of  force  on  the  perimeter  of  infinitely  long  moving  objects  calculated  from  FEM  simulation 
(—>•),  compared  to  predictions  from  RFT  obtained  by  simulation  (—>>).  Net  resistive  force  shown  at  the  center  of  the 
object. 


is  confined  due  to  no  slip  condition. 

Exploration  of  superposition  in  3D  -  Equation  1  can  be  extended  naturally  to  include  the  out- 
of-plane  dimension,  backed  by  experiments  [2]  that  have  verified  that  surfaces  whose  shapes  vary  in 
the  y  direction  also  maintain  the  superposition  principle.  To  assess  the  RFT  superposition  principle 
in  3D,  we  employ  our  continuum  model  to  study  a  sequence  of  horizontal  motions  of  a  buried  V- 
shaped  intruder.  There  are  natural  limitations  on  the  superposition  principle  regarding  the  curvature, 
angularity  or  orientation  of  the  intruding  object.  For  example,  if  the  V  angle  is  too  close  to  zero, 
i.e.  the  two  wings  of  the  V  nearly  coalesce,  RFT  will  be  wrong  by  a  factor  of  two  since  one  wing  is 
completely  blocked  by  the  other.  To  study  RFT  in  its  plausible  regime,  we  hence  limit  our  attention 
to  an  obtuse  V  geometry,  with  angle  fixed  at  135°  that  ought  to  be  within  the  limits  of  the  theory. 
The  intruder  is  located  at  depth  of  0.75m  measured  from  the  center  of  its  spine  and  each  wing  is  a 
0.2  x  0.2m  square. 

We  consider  the  case  where  the  intruder  is  moving  as  a  one  piece  object  and  we  record  its  forces 
in  three  dimensions.  Then  we  run  two  new  simulations,  one  for  each  wing  of  the  intruder,  where  the 
wing  is  simulated  as  a  separate  single  plate  maintaining  its  same  exact  same  alignment,  positioning, 
and  motion  that  it  had  when  it  partook  in  the  group  V  motion.  We  consider  both  vertically  and 
horizontally  aligned  intruder  cases,  varying  the  orientation  angle  <f>  denoting  a  pitching  angle  (in 
vertical  case)  or  yawing  angle  (in  horizontal  case),  see  Fig.  4.  Both  cases  engage  in  considerable  drag 
force  in  the  direction  of  motion  as  expected.  The  vertically  aligned  intruders  experiences  no  lateral 
force  due  to  symmetry,  with  lift  forces  that  switch  sign  due  to  the  plowing.  The  horizontal  intruder 
experiences  somewhat  steady  lift  force,  but  the  lateral  forces  vary  considerably  for  similar  lift  of  the 
vertical  intruder. 

We  evaluate  RFT  by  comparing  the  forces  obtained  directly  by  FEM  of  the  V  intruder  (Fig. 
4  solid  lines)  with  the  forces  obtained  by  summation  of  contributions  from  individually  moving 
component  plates  (Fig.  4  dashed  lines).  This  comparison  shows  that  the  theory  generally  works 


very  well  in  3D.  We  find  the  agreement  is  more  pronounced  when  a  particular  wing  is  positioned 
such  that  there  is  less  interaction  between  the  flow-fields  stemming  from  each  plate.  The  agreement 
is  roughly  in  the  same  range  as  the  deviations  observed  in  past  granular  RFT  experiments  [3,  14]. 
For  certain  orientations,  the  disagreement  is  due  to  a  ‘shadowing’  effect  caused  by  the  leading  wing 
on  the  trailing  wing.  For  instance  at  0°  a  good  agreement  between  group  and  individual  forces  is 
observed  for  the  vertically  aligned  intruder.  One  can  expect  the  flow  fields  of  both  wings  to  point 
in  the  same  direction  (See  Fig.  1)  and  little  effect  should  be  expected  from  the  cross-interaction 
of  the  fields.  But  once  the  same  object  is  rotated  by  90°,  these  fields  would  interact  and  therefore 
the  disagreement  in  forces  increases  in  all  components.  For  this  3D  study  there  are  no  analogous 
experimental  tests  in  the  literature  to  verify  the  FEM  findings.  But  the  results  shown  in  Fig.  4 
suggest  that  for  3D  objects  that  do  not  contain  sharp  acute  angles  and  their  geometry  does  not 
shadow  parts  of  their  own  perimeter,  reasonable  results  should  be  expected  if  RFT  is  used. 


9J 


0.3 

0.25 

0.2 

0.15 


—  0.1 

CD 

O 

*  0.05 


0 

-0.05 

-0.1 

-0.15 

-0.2 


x\  V* 

x  x/ 
V 


/' 

/  Fx 


/ 


/ 


Y'  v'\  \ 

^A-A-A-A-A-A-A^A -A-^ 

\  7  \ 

A  / 

\v  y 

V/ 


90 


180 

cp° 


270 


360 


Figure  4:  The  drag  force,  Fx,  acting  on  a  rightward  moving  submerged  ’V’  intruder  (x),  and  the  components  of  lift, 
Fz  (v  ),  and  lateral  force,  Fy  (A  ),  compared  with  superposition  values  (dashed)  at  various  orientation  angles,  ip,  for 
horizontal  (left)  and  vertical  (right)  intruder  alignments. 


Mathematical  explanation  -  The  results  thus  far  have  demonstrated  numerically  that  the  hy¬ 
potheses  of  RFT  emerge  relatively  strongly  from  the  continuum  equations  of  frictional  plasticity. 
To  provide  more  of  an  explanation  as  to  why  the  continuum  equations  replicate  RFT,  we  consider 
the  following  example.  As  delineated  earlier,  exact  solutions  to  the  continuum  plasticity  system  are 
highly  nontrivial  to  obtain,  but  the  example  we  highlight  below  can  be  proven  using  dimensional 
analysis  to  collapse  exactly  to  RFT. 

Consider  a  large,  quasi-2D  (of  arbitrary  out-of-plane  width  D ),  semi-infinite  bed  of  our  continuum 
plastic  media.  Suppose  a  straight  intruder  is  inserted  into  the  media  at  an  angle  /?  from  the  horizontal, 
such  that  L  of  its  length  is  submerged.  It  is  then  translated  at  a  speed  v  in  a  direction  angled  7 
from  horizontal,  producing  an  assumed  quasi-static,  plane-strain  motion  of  the  material.  At  the 
moment  the  motion  begins,  the  total  resistive  force  F  on  the  intruder  can  be  calculated.  If  the 


calculation  is  repeated  for  a  slightly  longer  intruder,  of  submerged  length  L  +  dZ,  we  can  ask  how 
much  extra  resistive  force  is  there,  or  equivalently,  what  is  dF /dll  The  continuum  model  inputs  and 
boundary  conditions  imply  that  the  answer  can  only  depend  on  /?,  7,  L,  the  granular  areal  weight 
density  pgD ,  the  internal  friction  /ic,  and  perhaps  a  wall  friction  coefficient  pw  between  the  intruder 
and  the  material.  Let  the  units  of  length,  time,  and  mass  be  [Length]  =  L,  [Time]  =  L/v  and 
[Mass]  =  pgD  *  L  *  (L/v)2  respectively.  The  dimensionless  groups  are  then  (dF  /  dl)  /  (pgD  *  L2),  /?, 
7,  pC:  and  pw.  It  therefore  follows  that 

dF 

—  =  pgDL  -0(/3, 7,  pc,  fj,w)  =  \z\D  ^(/3, 7,  pc,  gw)  (3) 

where  the  right  equality  uses  \z\  —  Lsin/3.  By  integrating,  we  fold  Eq.3  gives  the  general  form  of 
the  RFT  (see  Eq.l)  where  the  components  of  the  function  Vi/  are  precisely  the  lift  and  drag  functions 
7)  and  ctx(ft, 7),  which  depend  on  the  material  properties  and  intruder  roughness  through  pc 
and  pbw. 

While  instructive,  the  above  shows  only  that  RFT  is  exact  in  the  specific  case  described  and  does 
not  explain  its  general  effectiveness  for  curved,  submerged  surfaces.  In  more  general  geometries,  we 
hypothesize  that  the  reason  for  continued  RFT-type  superposition  may  be  linked  to  the  fact  that 
the  equations  of  frictional  plasticity  for  2D,  quasi-static,  rigid-plastic  problems  become  a  hyperbolic 
system  in  space.  The  stress  characteristics  extend  from  boundaries  along  lines  saturating  \r\/P  <  /ic, 
that  form  the  ‘slip-lines’  [15].  This  produces  domains  of  dependence  in  the  material  in  a  sense  that 
stresses  in  certain  zones  can  be  attributed  to  the  traction  on  a  specific  part  of  the  surface  of  the 
intruder.  This  is  in  sharp  contrast  to  viscous  fluids,  for  example,  in  which  the  equations  are  elliptic. 
Experimentally,  superposition  tends  to  work  better  in  granular  media  than  it  does  for  viscous  fluids 
[16],  which  is  particularly  ironic  given  that  RFT  was  initially  developed  for  viscous  drag  [17].  Using 
an  analogous  dimensional  analysis  as  above  but  this  time  assuming  a  viscous  fluid  rheology  leads  to 
dF/dl  =  |z|  D  Vi /(/?,  7,  gw,vz/ is)  which  confirms  that  the  form  of  the  force  function  does  not  reduce  to 
the  superposition  form  of  Eq.3  since  viscosity  brings  in  an  additional  velocity  and  depth  dependence 
through  a  dimensionless  group  vz/is  for  is  the  kinematic  viscosity.  Further  details  of  the  viscous  vs 
granular  analysis  can  be  found  in  the  Supplemental  Materials. 

We  have  shown  that  the  basis  of  the  granular  superposition  rule  and  the  other  hypotheses  of 
RFT,  arise  in  full  from  the  most  basic  law  of  continuum  frictional  plasticity.  We  have  verified  this 
point  directly  in  plane-strain  and  full  3D  problems.  Moreover,  the  plasticity  system,  when  ht  to  the 
repose  angle  of  the  physical  media  of  [2] ,  quantitatively  replicates  the  empirical  input  data  that  was 
previously  used  for  RFT.  We  have  identified  how  the  exact  form  of  RFT  arises  in  a  model  problem, 
based  solely  on  dimensional  analysis.  In  so  doing,  we  have  also  explained  why  RFT  works  better  in 
granular  matter  than  in  viscous  fluids.  This  new  fundamental  understanding  paves  the  path  forward 
to  better  understand  and  derive  RFT’s  for  other  flowable  materials,  or  apply  the  method  in  systems 
that  might  be  difficult  to  approach  otherwise  experimentally  such  as  the  modeling  of  terradynamics 
in  micro-gravity.  For  example,  using  a  similar  analysis  to  the  one  above,  certain  cohesive  media  may 
also  possess  an  RFT  describing  forces  on  moving  locomotors;  see  Supplemental  Materials  4.2.  By 
rooting  RFT  in  the  fundamental  plasticity  equations  we  now  have  a  clear  mechanical  framework  to 
understand  the  various  parameters  in  RFT  as  well  as  to  discern  the  circumstances  that  diminish  the 
accuracy  of  the  model. 
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The  goal  of  this  supplemental  material  is  to  provide  details  about  the  modeling  technique  and 
calculations  presented  in  the  main  text.  We  describe  different  components  of  the  finite  element 
approach  and  the  corresponding  parameters  followed  by  a  more  detailed  discussion  on  the  dimensional 
analysis. 

1  Numerical  solution  procedure 

1.1  Elasticity-augmented  constitutive  relation 

As  outlined  in  the  main  paper,  we  describe  the  deformation  behavior  of  the  granular  material  by  using 
a  non-hardening  Drucker-Prager  (DP)  yield  criterion  (constant  /ic).  When  implementing  the  model, 
we  assume  a  small  portion  of  the  deformation  is  elastic,  which  closes  the  system  mathematically  in 
regions  of  non-flowing  material,  and  provides  a  natural  route  to  implementing  the  pressure  constraints 
described.  The  model  then  takes  the  form  of  a  simple  hypo-plastic  formulation.  At  any  time,  the 
deformation  gradient  F{j  is  used  to  construct  the  velocity  gradient  dvi/dxj  —  which  is  further 

divided  into  a  symmetric  stretching  component  Dij  and  anti-symmetric  spin  component  W{j.  The 
dot  represents  the  Lagrangian  time  derivative.  The  stretching  is  further  decomposed  into  elastic  and 
plastic  parts  =  D\-  +  IF-.  Denoting  the  elastic  shear  and  bulk  moduli  of  the  granular  material 
by  G  and  K  respectively,  the  constitutive  relationship  for  the  stress  evolution  of  dense  material  (i.e. 
p  >  Pc)  is  established  by  a  Jaumann  objective  rate  form  as  follows: 

I °ij  =0,  if  P  <  Pc  ^ 

\&ij  -Wik(Tkj  +  VikWkj  =  KDlk5ij  +  2G(D?j  -  Dlkdij/3),  iip>pc 

where  <Jij  is  the  Cauchy  stress  as  described  in  the  main  paper  which  is  divided  into  a  hydrostatic 
pressure  part  a  —  cr^h/ 3  =  —  P  and  a  deviatoric  part  a-  =  <Jij  —  aftij.  This  relationship  warrants 
that  the  material  is  stress-free  if  density  falls  below  a  certain  critical  density  pc  which  is  indicative 
of  a  disconnected  or  ‘open  state’  of  the  material.  Further  discussion  about  this  condition  is  followed 
in  section  1.3  of  this  supplement. 

1.2  Update  step 

The  numerical  scheme  is  a  variant  of  the  radial  return  algorithm  [1]  .  It  starts  each  step  by  assuming 
a  purely  elastic  step,  i.e.  D^-  —  0,  under  Equation  1,  which  updates  the  stress  to  a  “trial  stress” 

(crF).  If  the  trial  stress  results  in  an  equivalent  shear  stress  ftr  =  \J Gtj' that  ^ess  than  PcPtr , 
it  is  accepted  as  the  updated  stress.  If  not,  it  is  then  adjusted  as  shown  in  Figure  1.  In  this  figure, 
a\-  is  the  stress  state  at  the  beginning  of  the  step.  Since  Ptr  —  Pf+At  due  to  the  isochoric  plastic  flow 


assumption,  the  effective  shear  stress  f  should  reduce  following  a  constant  pressure  route  to  reside 
on  the  yield  surface  and  cr^At,  the  stress  state  at  the  end  of  increment,  is  updated  accordingly.  This 
essentially  represents  usage  of  a  tangent  modulus  to  return  the  trial  stress  state  ajj  back  to  the  yield 
surface  at  the  end  of  the  increment  [1]. 


Figure  1:  Projection  of  stresses  to  the  yield  surface. 


4+At  =  (ncPt+At/ftr)  -  Pt+A%  (2) 

Note  that  the  values  of  P  and  f  should  be  updated  after  this  projection.  Adjustments  should  also 
be  made  if  the  material  is  stretched  beyond  a  critical  density  as  described  below. 

1.3  Simultaneous  dense  and  disconnected  matter 

Following  [2],  the  opening  rule  is  implemented  by  tracking  the  Jacobian  (det(i^j))  at  each  integration 
point.  If  the  initial  local  density  at  a  material  point  is  denoted  by  po5  then  its  density  in  deformed 
state  is  p  —  po/det(F^).  If  the  Jacobian  exceeds  unity,  the  local  density  at  the  deformed  stage  is 
reduced.  Once  the  deformed  density  is  smaller  than  critical  density  pc  the  material  point  is  flagged 
to  behave  locally  as  open  state  i.e.  stress  free  deformation.  In  all  our  simulations,  we  begin  by 
assuming  po  =  Pa  he.  a  granular  assembly  barely  in  contact,  and  turn  gravity  on  in  a  following  step, 
which  causes  the  material  to  compress  elastically  a  small  amount. 

In  the  case  of  an  flat  plate  traveling  within  granular  material,  one  would  expect  that  behind 
the  intruder,  the  free  space  caused  by  movement  of  the  object  would  cause  open  states  in  its  wake. 
Consequently  the  open  material  above  the  region  should  fall  down  and  fill  the  cavities  caused  by  the 
intruder  movement  as  shown  in  density  plot  in  2.  The  open  states  in  the  back  of  the  intruder  are  more 
or  less  persistent  after  reaching  steady-state  deformation  while  in  the  front  there  are  isolated  areas 
that  show  a  limited  amount  of  reduction  in  density.  As  expected,  the  figure  shows  these  locations 
form  in  the  trailing  path  of  the  intruder. 


Figure  2:  Plot  of  the  local  density  in  the  granular  bed  in  the  wake  of  a  flat  vertical  intruder  moving  rightward  using 
FEM.  Densities  below  pc  =  2000  kg/m3  denote  disconnected  material  in  an  ’open  state’. 


2  Dimensional  analysis 

Here  we  extend  and  clarify  the  dimensional  analysis  examples  presented  in  the  main  paper  with  a 
complete  description  to  study  the  form  of  resistive  forces  in  granular  materials  and  to  contrast  it 
with  viscous  fluids. 


2.1  Granular  material 

Assume  a  perfectly  plastic  bed  of  grains  (AT,  G  — >>  oo)  and  a  straight  intruder  of  out-of-plane  width 
D  and  length  L.  Let  the  friction  coefficient  for  the  granular  media  and  the  wall  friction  coefficient 
for  the  intruder  be  (ic  and  /iw  respectively,  and  the  weight  density  of  the  material  be  pg.  Assume 
quasi-static  motion  of  the  intruder. 


2.1.1  Intruder  in  the  vertical  plane 

Consider  a  cross-sectional  slice  of  the  above  mentioned  box  as  shown  in  figure  3  in  plane  strain 
conditions.  The  intruder  travels  at  speed  v  with  intrusion  angle  of  7,  pushing  the  material  in  front. 
The  intruder  is  oriented  at  angle  /3  with  its  length  of  L  being  in  contact  with  the  granular  material. 
We  assume  plasticity  is  reached  instantaneously.  If  the  force  is  denoted  as  F,  then  the  rate  of  change 
in  force  with  respect  to  change  in  length,  denoted  by  dF/dZ,  represents  the  deviation  in  the  force  if 
the  intruder  were  extended  slightly  longer  by  length  of  dZ.  Based  on  the  assumptions  above,  dF/dZ 
can  only  depend  on  /?,  7,  L,  p,  pgD,  pc  and  pw.  The  units  of  length,  time  and  mass  are  represented 
as  [L\  =  L,  [T]  =  L/v,  [M]  =  pgD  *  L  *  L2 /v2  respectively  and  the  dimensionless  force  denoted  by 
tilde  is  obtained  as: 

F  F  F 

F  =  PPI7PF  =  =  ^Di?  (3) 

Therefore  it  follows  that  dF/dl  =  pgDL  7,  pc,  pw).  Substituting  \z\  =  L  sin  /?  ,  we  can  write: 
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Figure  3:  Arbitrarily  oriented  vertical  intruder  set  into  motion  in  an  arbitrary  direction.  Plane  strain  conditions. 


which  is  the  differential  form  of  the  RFT  equation;  pg  has  been  absorbed  in  to  the  definition  of  \l/ vert . 
Here  the  components  of  the  function  ^vert  precisely  define  ax(/3, 7)  and  07  (/3, 7)  of  Eq  1  in  the  main 
text. 

2.1.2  Horizontal  intruder 

Consider  the  same  configuration  described  above,  but  this  time  for  a  horizontally  aligned  intruder. 
This  time,  we  consider  a  horizontal  cross-section  at  depth  z  >  0  passing  through  the  intruder.  By 
making  the  simplifying  approximations  that  (i)  that  the  pressure  is  a  linear  function  of  depth,  and 
(ii)  flow  is  contained  within  the  cross-section,  we  can  identify  a  pseudo-cohesion  parameter  in  the 
equations,  defined  as  pcpgzD ,  which  emerges  as  a  constant  effective  material  strength  per  length 
within  the  cross-section.  Therefore  dF/dl  in  this  case  depends  on  9,lj,L,v,  pcpgzD,  pc,  and  pw.  We 
choose  dimensional  units  in  this  configuration  to  be  [L\  =  L,  [T]  —  L/v,  [M]  —  gcpgzD  *  L2/v2  and 
the  dimensionless  force  is: 


F  F  F 

[M][L]/[s]2  p>cpg  z  L  VcpgzL 

Therefore  it  follows  similarly  that: 


dF 

~dl 


=  ncpgzD  iphor(9,uj,pw,pc) 


dF 

~dl 


Z  |  D  ^ /ior(^  5^  1  pwi  pc) 


(5) 

(6) 


Figure  4:  Horizontal  intruder  -  plane  view  perpendicular  to  gravity. 

This  analogy  creates  a  correlation  between  the  effects  of  depth  in  cohesion-less  material  and 
arbitrary  cohesion  parameter  in  a  cohesive  material. 


2.2  Stokes  flow:  Intruder  in  a  vertical  plane 


For  viscous  incompressible  slow-regime  flow,  the  general  Navier-Stokes  equation  reduces  to  Stokes 
equation  as  —np(d2/dx2)vi  —  ( d/dxi)P  +  fi  —  0  with  condition  dvi/dxi  —  0  to  enforce  incompress¬ 
ibility.  Here,  v  is  velocity,  P  is  pressure,  v  is  kinematic  viscosity  and  /  is  body  force.  For  the  case 
of  vertical  intruder,  the  only  body  force  is  the  gravity  and  the  boundary  condition  and  properties  of 
the  material  imply  that  the  forces  depend  on  <a,  /?,  v,  L,  z,  pgD ,  v.  The  main  units  are  again  found  as 
[L\  —  L,  [T]  =  L/v,  [M]  —  pgD  *  L  *  ( L/v )2,  then  the  dimensionless  viscosity  is  found  by  v  =  v/vL 
and  the  dimensionless  force  is: 


F  F  F 

F  =  WmJW  =  =  pgV  <7) 

Consequently,  the  force  on  the  intruder  is  found  as: 

dF  v  dF  v 

—  =  pgL  ipvisc{f3, 7,Mw»  —f)  ~rr  =  \z\  ^ visc(P , 7 ,  ,  — )  (8) 

Cll  V  ij  cU  vz 

As  seen  above,  the  dimensionless  viscosity  group  carries  z-dependence  inside  the  function  \I &Visa- 

This  dependence  enables  the  correlation  to  deviate  from  a  purely  linear  correspondence  on  z  as  seen 
in  RFT  for  granular  material.  A  similar  result  occurs  for  the  case  where  the  object  is  horizontally 
oriented. 
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